Formation of a vortex lattice in a rotating BCS Fermi gas 
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We investigate theoretically the formation of a vortex lattice in a superfluid two-spin component 
Fermi gas in a rotating harmonic trap, in a BCS-type regime of condensed non-bosonic pairs. Our 
analytical solution of the superfluid hydrodynamic equations, both for the 2D BCS equation of 
state and for the 3D unitary quantum gas, predicts that the vortex free gas is subject to a dynamic 
instability for fast enough rotation. With a numerical solution of the full time dependent BCS 
equations in a 2D model, we confirm the existence of this dynamic instability and we show that it 
leads to the formation of a regular pattern of quantum vortices in the gas. 
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The field of trapped ultracold fermionic atomic gases 
is presently making rapid progress: thanks to the pos- 
sibility of controlling at will the strength of the s-wave 
interaction between two different spin components by the 
technique of the Feshbach resonance PI H| , it is possible 
to investigate the cross-over Q between the weakly inter- 
acting BCS regime (case of a small and negative scatter- 
ing length) and the Bose-Einstein condensation of dimers 
(case of small and positive scattering length), includ- 
ing the strongly interacting regime and even the unitary 
quantum gas (infinite scattering length) . The interaction 
energy of the gas was measured on both sides of the Fes- 
hbach resonance 0; on the side of the resonance with 
a positive scattering length, Bose-Einstein condensation 
of dimers was observed ; on the side of the resonance 
with a negative scattering length, a condensation of pairs 
was revealed in the strongly interacting regime by a fast 
ramping of the magnetic field across the Feshbach reso- 
nance [5(- Also, the presence of a gap in the excitation 
spectrum was observed @ , for an excitation consisting in 
transferring atoms to an initially empty atomic internal 
state, as initially suggested by pfl, revealing pairing. 

Are there evidences of superfluidity in these fermionic 
gases ? It was initially proposed [j| to reveal superfluidity 
by detecting an hydrodynamic behavior in the expansion 
of the gas after a switching-off of the trapping potential. 
Such an hydrodynamic behavior was indeed observed 
but it was then realized that this can occur not only in 
the superfluid phase, but also in the normal phase in the 
so-called hydrodynamic regime, that is when the mean 
free path of atoms is smaller than the size of the cloud, 
a condition easy to fulfill close to a Feshbach resonance. 
The general experimental trend is now to try to detect 
superfluidity via an hydrodynamic behavior that has no 
counterpart in the normal phase 9J. A natural candi- 
date to reveal superfluidity is therefore the detection of 
quantum vortex lattices in the rotating trapped Fermi 
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gases: the superfluid velocity field, defined as the gradi- 
ent of the phase of the order parameter, is irrotational 
everywhere, except on singularities corresponding to the 
vortex lines, so that a superfluid may respond to rotation 
by the formation of a vortex lattice jlfj ; on the contrary, 
a rotating hydrodynamic normal gas is expected to ac- 
quire the velocity field of solid-body rotation and should 
not exhibit a regular vortex lattice in steady state. 

Steady state properties of vortices in a rotating Fermi 
gas described by BCS theory have already been the sub- 
, ect of several studies, for a single vortex configuration 
111 and more recently for a vortex lattice configuration 
1Q| . In this paper, we study the issue of the time de- 
pendent formation of the lattice in a rotating Fermi gas, 
by solving the time dependent BCS equations. A cen- 
tral point of the paper is to identify possible nucleation 
mechanisms of the vortices that could show up in a real 
experiment. 

This problem was addressed a few years ago for rotat- 
ing Bose gases. The expected nucleation mechanism was 
the Landau mechanism, corresponding to the apparition 
of negative energy surface modes for the gas in the rotat- 
ing frame, for a rotation frequency above a minimal value; 
these negative energy modes can then be populated ther- 
mally, leading to the entrance of one or several vortices 
from the outside part of the trapped cloud 0, ^| . The 
first experimental observation of a vortex lattice in a ro- 
tating Bose-Einstein condensate revealed however a nu- 
cleation frequency different from the one of the thermal 
Landau mechanism and was suggested later on to 
be due to a dynamic instability of hydrodynamic nature 
triggered by the rotating harmonic trap | 1 5 [. which was 
then submitted to experimental tests |l6l Il7j . The dis- 
covered mechanism of dynamic instability was checked, 
by a numerical solution of the purely conservative time 
dependent Gross-Pitaevskii equation, to lead to turbu- 
lence 0| and to the formation of a vortex lattice [T^ |. 

We now transpose this problem to the case of a two 
spin component Fermi gas, initially at zero temperature 
and stirred by a rotating harmonic trapping potential of 
slowly increasing rotation speed, as described in section 
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[I] Does the hydrodynamic instability phenomenon oc- 
cur also in the fermionic case, and does it lead to the 
entrance of vortices in the gas and to the subsequent for- 
mation of a vortex lattice ? We first address this prob- 
lem analytically, in section^ by solving exactly the time 
dependent two-dimensional hydrodynamic equations and 
by performing a linear stability analysis: very similarly 
to the bosonic case, we find that a dynamic instability 
can occur above some minimal rotation speed. We also 
extend this conclusion to the 3D unitary quantum gas. 
Then we test this prediction by a numerical solution of 
the time dependent BCS equations on a two-dimensional 
lattice model, in section ITTT1 this confirms that the dy- 
namic instability can take place and leads to the entrance 
of vortices in the gas, which are then seen to arrange in 
a regular pattern at long evolution times. 



I. OUR MODEL 

We consider a gas of fermionic particles of mass to, 
with equally populated two spin states j and J,, trapped 
in a harmonic potential and initially at zero tempera- 
ture. The particles with opposite spin have a s-wave 
interaction with a negligible range interaction potential, 
characterized by the scattering length a^o, whereas the 
particles in the same spin state do not interact. 

We shall be concerned mainly by the limit of a 2D 
Fermi gas. In this case, the trapping potential is very 
strong along z axis so that the quantum of oscillation 
along z, that is Hu! z , where tu z is the oscillation frequency 
along z, is much larger than both the mean oscillation 
energy in the x — y plane and the interaction energy per 
particle, so that the gas is dynamically frozen along z in 
the ground state of the corresponding harmonic oscilla- 
tor. In this geometry, the two-body interaction can be 
characterized by the 2D scattering length a 2 D which was 
calculated as a function of the 3D scattering length in 
[20|. We recall that a 2 D is always strictly positive and 
the 2D two-body problem in free space exhibits a bound 
state, that is a dimer, of spatial radius a 2 o- For the 
2D gas to have universal many-body interaction proper- 
ties, characterized by a2D only, one requires that the spa- 
tial extension (h/mwz) 1 ^ 2 of the ground state of the har- 
monic oscillator along z is smaller than a 2 u HD> so that 
e.g. the dimer binding energy is smaller than fiw 2 . The 
weakly attractive Fermi gas limit corresponds in 2D to 
P a 2D ~~ * + 00 an d the condensation of preformed dimcrs 
to pa\ D — ► [23, where p is the 2D density of the gas. 

In the x — y plane, the zero temperature 2D gas 
is initially harmonically trapped in the non-rotating, 
anisotropic potential 



potential into rotation around z axis with an instanta- 
neous rotation frequency until it reaches a maximal 
value O to which it then remains equal. The question is 
to study the resulting evolution of the gas and predict 
the possible formation and subsequent crystallization of 
quantum vortices. 

We shall consider this question within the approxi- 
mate frame of the BCS theory, in a rather strongly inter- 
acting regime but closer to the weakly interacting BCS 
limit than to the BEC limit, which is most relevant for 
the present 3D experimental investigations: the chemi- 
cal potential p of the 2D gas is supposed to be positive, 
excluding the regime of Bose-Einstein condensation of 
the dimers, and the parameter kpci 2 D, where the Fermi 
momentum is defined as h 2 kp/2m = p, is larger than 
unity but not extremely larger than unity: we shall take 
= 4 in the numerical simulations. In this rela- 
tively strongly interacting regime, we of course do not 
expect the BCS theory to be 100% quantitative. 

In the hydrodynamic approach to come, one simply 
needs the equation of state of the gas, that is the expres- 
sion of the chemical potential pa of a spatially uniform 
zero temperature gas as a function of the total density 
P = Pi + Pi = 2/0f and of the scattering length. In 2D, 
this equation of state was calculated with the BCS ap- 
proach in p2| : 
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where Eg is the binding energy of the dimer in free space, 
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and 7 = 0.57721 ... is Euler's constant. Similarly, the 
gap for the zero temperature homogeneous BCS gas is 
related to the density by |22j 



Ao[p] = E a 



2nh 2 p 
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We shall also consider analytically the 3D unitary quan- 
tum gas (CL3D — 00) where the equation of state is known 
to be exactly of the form po[p] oc h 2 p 2 / 3 /m. 

In the numerical solution of the 2D time dependent 
BCS equations to come, one needs an explicit micro- 
scopic model. We have chosen a square lattice model 
with an on-site interaction between opposite spin parti- 
cles corresponding to a coupling constant go so that the 
second quantized grand canonical Hamiltonian reads at 
the initial time 



U(r) = \mu 2 [(1 



e)x 2 + (l + e)y 2 
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where r = (x,y) and e > measures the anisotropy of the 
trapping potential. Then one gradually sets the trapping 
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where I is the grid spacing. In the numerics a quantiza- 
tion volume is introduced, in the form of a square box 
of size L with periodic boundary conditions, L being an 
integer multiple of I. The sum over r then runs over the 
(L/l) 2 points of the lattice. A plane wave on the lattice 
has wavevector components k x and k y having a meaning 
modulo 2ir/l so that the wavevector k is restricted to the 
first Brillouin zone D = [—ir/l,Tr/l[ 2 . The operator cu,a 
annihilates a particle of wavevector k and spin state a —\ 
or I, and obeys the usual fermionic anticommutation re- 
lations, such as 



{Ck,, 



ji j <5k,k"5<T,a 
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The discrete field operator is proportional to the 

annihilation operator of a particle at the lattice node r 
in the spin state a in such a way that it obeys the anti- 
commutation relations 



{Mr),^l(r')} = r 2 S ry S a 



(7) 



The coupling constant go is adjusted so that the 2D scat- 
tering len gth of two particles on the infinite lattice is ex- 
actly a 2 i)l2i|23: 
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where G = 0.91596 ... is Catalan's constant. In the limit 
a- 2 D ~ * +oo, for a fixed density p and a fixed 'range' 
I of the interaction potential, one finds go — > CP: we 
recover the fact that the limit kp^D ^> 1 corresponds to 
a weakly attractive Fermi gas. 

At later times, the Hamiltonian is written in the frame 
rotating at frequency to eliminate the time depen- 
dence of the trapping potential; this adds an extra term 
to the Hamiltonian. 
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where the matrix L z on the lattice represents the angular 
momentum operator along z, xp y — yp x - The square box 
defining the periodic boundary conditions is supposed to 
be fixed in the rotating frame, so that it rotates in the 
lab frame: this may be useful in practice to ensure that 
truncation effects due to the finite size of this box in the 
numerics do not arrest the rotation of the gas. 

This lattice model is expected to reproduce a contin- 
uous model with harmonic trapping and zero range in- 
teraction potential in the limit of an infinite quantization 
volume (L spatial radius of the cloud) and in the limit 
of a vanishing grid size I — > (I <C a,2D, kp ). In this limit 
go is negative, leading to an attractive interaction, so that 
pairing can take place in the lattice model. In this limit, 
we have checked that BCS theory for the lattice model 
gives the same equation of state as Eq.@ 2j|. 



II. SOLUTION TO THE SUPERFLUID 
HYDRODYNAMIC EQUATIONS 

In the hydrodynamic theory of a pure superfluid with 
no vortex, one introduces two fields, the total spatial den- 
sity of the gas, p(r, t), and the phase of the so-called order 
parameter, 2 S(r,t)/h. In the BCS theory for the lattice 
model, the order parameter is simply 

A(r,t) ee -g (^(r,t)^(r,t)) = \A\e 2lS / h (10) 

which has a finite limit when I — > 0. The superfluid 
velocity field in the lab frame is then defined as 

v = M^. (ii) 

TO 

In the rotating frame, the hydrodynamic equations 
read 



d t p = — div [p (v — Q(t) 



x r 



-d t S = -mv 2 + U(r)+ Mo[p(r,t)] 
— fi — m(fl(t) x r) ■ v 



(12) 
(13) 



where S~2(t) = f2(<)z and z is the unit vector along the ro- 
tation axis z. The first equation is simply the continuity 
equation in the rotating frame, including the fact that the 
velocity field in the rotating frame differs from the one in 
the lab frame by the solid body rotational term. When 
one takes the gradient of the second equation, one recov- 
ers Euler's equation for a superfluid. These superfluid 
equations are expected to be correct for a slowly varying 
density and phase, both in space (as compared to the size 
of a BCS pair) and in time (as compared to ft/|A|). For 
a harmonically trapped system with a quantum of oscil- 
lation fku, the slow spatial variation condition requires a 
gap parameter |A| ^> ftw: in the present paper, consid- 
ering the rather strongly interacting regime 1 < kpa-iD-, 
the gap is of the order of the Fermi energy, which is much 
larger than fku, so that there is slow spatial variation as 
long as no vortex enters the cloud. The gap is then much 
larger than h over the ramping time of the trap rotation, 
so that the expected condition of slow time variation is 
also satisfied. In the appendix^] we present a simple but 
systematic derivation of these superfluid hydrodynamic 
equations starting from the time dependent BCS theory 
and using a semi-classical expansion. Surprisingly, for 
the case of slow ramping times and rather fast rotations 
considered in this paper, with f2 of the order of uj, our 
simple derivation requires an extra validity condition, in 
general more stringent than |A| 3> hco: the quantum of 
oscillation ftw should be smaller than |A| 2 //i, a condition 
also satisfied in our simulations. 

We shall assume here that the rotation frequency is 
ramped up very slowly so that the density and the phase 
adiabatically follow a sequence of vortex free stationary 
states. The strategy then closel y fo llows the one already 
developed in the bosonic case [l5|]: one solves analyti- 
cally the corresponding stationary hydrodynamic equa- 
tions, then one performs a linear stability analysis of the 



4 



stationary solution. The apparition of a dynamic insta- 
bility suggests that the system may evolve far away from 
the stationary branch; that this dynamic instability re- 
sults in the entrance of vortices will be checked by the 
numerical simulations of section ITTT1 

In the stationary regime, for a fixed rotation frequency 
fl one sets d t p = in Eq.JUJl and -d t S = in Ea. (jT3|) 
[26| . We first consider the 2D case and we replace (j,q 
by the equation of state Eq.Q: apart from an additive 
constant, (j,q is proportional to the density, as was the 
case for the weakly interacting condensate of bosons , 
so that the calculations for the superfluid fermions are 
formally the same, if one replaces the coupling constant 
g of the bosons by TrH 2 /m. Since the properties of the 
bosons do not depend on the value of g up to a scaling 
on the density 15], the results for the bosons can be 
directly transposed. Following [27]], we take the ansatz 
for the phase: 



S(r) = mujf3xy 



(14) 



which is applicable for a harmonic trapping potential 
U. When inserted in Eci. (|13[) . this leads to an inverted 
parabola for the density profile, resulting in an elliptic 
boundary for the density of the cloud. Upon insertion of 
the density profile in the continuity equation, one recov- 
ers the cubic equation of p?! ]: 



/? 3 



9? 



9. 



1-2— ) 0-e- =0. 



(15) 



This equation has one real root for 91 below some e depen- 
dent bifurcation value Jlbif(e), and has three real roots 
for 91 > 9 bi f(e). In the considered stirring procedure, 
the system starts with (3 — and follows adiabatically 
the so-called upper branch of solution, corresponding to 
increasing values of (3. In figure^] we have plotted (5 as a 
function of 91 /to on this branch, for the value of the asym- 
metry parameter in the simulations of the next section, 
e = 0.1. When (3 takes appreciable values, the cloud sig- 
nificantly deforms itself in real space, becoming broader 
along x axis than along y axis, even for an arbitrarily 
weak trap anisotropy e. 

From the studies of the bosonic case 0] it is known 
that the significantly deformed clouds can become dy- 
namically unstable. We recall briefly the calculation pro- 
cedure: one introduces initially arbitrarily small devia- 
tions bp and SS of the density and the phase from their 
stationary values; one then linearizes the hydrodynamic 
equations Eci. ((T2*|) and Eci. ((T3j> to get 



DSp 
~Dt 
DSS 
Dt 



-div I p 



grad SS 



dp 



(16) 
(17) 



where D /Dt = dt + (v — 91 x r) • grad and where we used 
the fact that the Laplacian of S(r) cx xy vanishes. One 
then calculates the eigenmodes of the linearized equa- 
tions, setting dt — > —iv where v is the eigenfrequency of 



CO. 




0.4 0.6 



FIG. 1: The upper branch of solution for the phase parameter 
/3 of the hydrodynamic approach for a stationary vortex free 
BCS state in the rotating frame, as a function of the rotation 
frequency. Solid line: the trap anisotropy is e = 0.1. Dashed 
line: e = 0. 



the mode. As an ansatz for Sp(r) and 8S(r), one takes 
polynomials of arbitrary total degree n in the variables x 
and y. One can indeed check that the subspace of polyno- 
mials of degree < n is stable, since the stationary values 
p and S are quadratic functions of x and y. This turns 
the linearized partial differential equations into a finite 
size linear system whose eigenvalues can be calculated 
numerically. Complex eigenfrequencies, when obtained, 
lead to a non-zero Lyapunov exponent A = Im v, which 
reveals a dynamical instability when A > 0. 

In figure |21 we plot the stability diagram of the upper 
branch stationary solution in the plane (91, e), for various 
total degrees n of the polynomial ansatz. Each degree 
contributes to this diagram in the form of a crescent, 
touching the horizontal axis (e — 0) with a broad basis 
on the right side and a very narrow tongue on the left side 
pSjj . For the low value e = 0.1 considered in the numer- 
ical simulations of this paper, the Lyapunov exponents 
in the tongues are rather small, so that significant insta- 
bility exponents are found only in the broad bases: for 
increasing 91, the first encountered significant instability 
corresponds to a degree n = 3: for e = 0, the correspond- 
ing minimal value of 9L/u is ( 183 +^ VM ) 1 / 2 = 0.79667 . . . 
[30j. This is apparent in figure where we plot the Lya- 
punov exponent as a function of 91 /lu for various degrees 
n and for e = 0.1. 

Extension to the unitary quantum gas in 3D: In practice, 
the experiments are mainly performed in 3D, so that we 
generalize the previous hydrodynamic calculation to a 
3D case where the exact equation of state is known: the 
so-called unitary regime, where the 3D s-wave scattering 
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FIG. 2: For the upper branch of solution for the phase pa- 
rameter, in 2D: Dark areas: instability domain in the Q — e 
plane for degrees n equal to 3, 4 and 5 (crescents from bottom 
to top). There is no dynamical instability for n < 2. Solid 
line: border Q, 2 — (1 — e)ui of the branch existence domain. 



length between opposite spin fermions is infinite. Be- 
cause of the universality of the unitary quantum gas, the 
equation of state of the gas is indeed a power law 



p [p] = Ap< 



(18) 



where the exponent 7 = 2/3 and where the factor A 
is proportional to h 2 /m, with a proportionality constant 
recently calculated with fixed node Monte Carlo methods 
[3lj, l32j and measured in recent experiments by Grimm 
|33| and by Salomon Q . 

For such a non-linear equation of state, one seems to 
have lost the underlying structure of the hydrodynamic 
equations allowing a quadratic ansatz for p and S, and 
a polynomial ansatz for 5p and SS. Fortunately, this 
structure can be restored by using as a new variable 
R(r,t) = p 7 (r, t). One then gets effective hydrodynamic 
equations with a linear equation of state: 

d t R = -jRdivv- (v-fl(t) x r) gradi? (19) 
-8 t S = ^mv 2 + U 3D (r) + AR(r) 

-p - m(0(i) x r) • v, (20) 

where the 3D trapping potential is 

U 3D (r) = X -mu? [(1 - e)x 2 + (1 + e)y 2 } + ^rmj 2 z z 2 . 

One then recycles the previous approach, with the usual 
quadratic ansatz for the steady state values of R and S. 



FIG. 3: For the upper branch of solution for the phase param- 
eter in 2D: Lyapunov exponent of the dynamic instability for 
degrees n from 3 to 7, as a function of the rotation frequency. 
The trap anisotropy is e = 0.1. 



In particular the same cubic equation for (3 as in Eq. H15(l 
is obtained. Linearizing the effective hydrodynamic equa- 
tions around the steady state, one gets 



DSR 

Dt 
DSS 

Dt 



- 7 i?^^-— grad(5S -gradi? (22) 
m m 



= -A5R, 



(23) 



where we used the fact that S has a vanishing Laplacian. 
This system of partial different equations can be solved 
by a polynomial ansatz for SS and SR. This generalizes 
to the rotating case the ansatz of |34| . 

In figure 21 we have plotted the stability diagram of 
the upper branch stationary solution in the plane (Q, e) 
for the 3D unitary quantum gas, for a trapping potential 
with lu z — OAuj. The 3D nature of the problem makes the 
structure of the instability domain more involved that in 
2D. This also appears in figure |SJ giving the Lyapunov 
exponents as a function of Q for a fixed trap anisotropy in 
the x — y plane, e = 0.022. In the limit of a cigar shaped 
potential, oj z <C lu, the structure is on the contrary close 
to the 2D one, as some of the eigenmodes for SR and SS 
almost factorize in a function of x, y and a function of z. 



III. NUMERICAL SOLUTION OF THE 2D 
TIME DEPENDENT BCS EQUATIONS 

We recall briefly the BCS equations for our two- 
component lattice model, in the case of equal popula- 
tions of the two spin states. In the non-rotating case, 
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FIG. 4: Case of the 3D unitary quantum gas with uo z = QAuj, 
for the upper branch of solution for the phase parameter: 
Dark areas: instability domain in the Q. — e plane for degrees 
(a) n = 3, (b) n = 4, (c) n = 5 and (d) n = 6. There is no 
dynamical instability for n < 2. 



the many-body ground state of the Hamiltonian is ap- 
proximated variationally in the zero temperature BCS 
theory by a so-called quasiparticle vacuum |35| . that is 
the vacuum state of annihilation operators of elementary 
excitations, 6 SjCT (where a =\ or J,). By energy minimiza- 
tion, one finds that the b s ^ a are such that 



V T (r) - J2[b sA u s (v)-bl d v* s (v) 

S 

V^i(r) = J2[b Sll u s (r) + bl 7 v* 3 (r) 



(24) 
(25) 



where the u's and v's are all the eigenvectors of the fol- 
lowing Hermitian system with positive energies E s > 0: 



h 
A* 



A 

-hi 



and normalized so that 



[| Us (r)| 2 + k(r)| 5 



(26) 



(27) 



In the eigensystem, A is the position dependent gap 
parameter defined in Ea. (|10f> and the matrix ho repre- 
sents on the lattice the single particle kinetic energy plus 
chemical potential plus harmonic potential energy terms. 
When the modal decompositions Ecis. 124I25H are inserted 
in Eq. H10[l . one gets 



A(r) = -g ^2u s (r)v*(r) 



(28) 



FIG. 5: Case of the 3D unitary quantum gas with uj z = 0.4w, 
for the upper branch of solution for the phase parameter: 
Maximal Lyapunov exponent of the dynamic instability for 
degrees n from 3 to 6, as a function of the rotation frequency. 
The trap anisotropy is e = 0.022. 



The density profile of the gas is given by 
p(r)=2(^(r)^ T (r)} = 2^K 



(29) 



These equations actually belong to the zero temperature 
Hartree-Fock-Bogoliubov formalism for fermions and are 
derived in §7. 4b of [3f|. Note that we have omitted the 
Hartree-Fock mean field term [3|j. 

To solve numerically the 2D self-consistent stationary 
BCS equations, we have used the following iterative al- 
gorithm: one starts with an initial guess for the position 
dependence of the gap parameter (we used the local den- 
sity approximation, taking advantage of the fact that the 
equation of state Eq.@ and the value of the gap Eq.QJ 
within BCS theory are known analytically in 2D), then 
one calculates the u's and v's by diagonalization of the 
Hermitian matrix in Ea. H26l) . one calculates the corre- 
sponding A(r) using Ea. (|28|l . and one iterates until con- 
vergence. 

Once the stationary BCS state is calculated, one moves 
to the solution of the 2D time dependent BCS equa- 
tions, to calculate the dynamics in the rotating trap. 
What we call here time dependent BCS theory is the 
time-dependent Hartree-Fock-Bogoliubov formalism for 
fermions, in the form of a variational calculation with 
a time dependent quasiparticle vacuum \(f>(t)), as de- 
tailed in §9.5 of [35j]. At time t, the modal expansions 
Eas. l|24l25|l still hold for ipi(r) and V'i( r ) I except that 
the operators (where a —\ or j) and the mode func- 
tions are now time dependent. The variational state vec- 
tor \cf>(t)) is the vacuum of all the operators 6 Sj(T (i). The 
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mode functions evolve according to 

«* (::)=(?. _«)(::) <3o » 

where ho now includes the rotational term — Q(t)L z in 
addition to the kinetic energy, the chemical potential and 
the trapping potential. The gap function A is still given 
by Eq.lEU and is now time dependent as the mode func- 
tions are. Note that Ea. l|3U|l corresponds to the first of 
the equations (9.63b) in §9.5 of |3^|, up to a global com- 
plex conjugation. To be complete, we give the expression 
of the time dependent quasiparticlc annihilation opera- 
tors: 

6 a , t (t) - Z 2 ^<(r,i)V> T (r) + <(r,t)^(r) (31) 

r 

b, ti (t) = Z 2 ^<(r,t)^(r)-<(r,t)^](r). (32) 

r 

We also recall that this time-dependent formalism con- 
tains not only pair-breaking excitations, but also implic- 
itly collective modes of the gas, as can be shown by a 
linearization of these equations around a steady-state so- 
lution, see §10.2 in [35|, and as also shown by the fact 
that hydrodynamic equations may be derived from them 
as done in the Appendix^] The numerical simulations 
to come therefore include excitations of these collective 
modes, when the numerical solution deviates from a sta- 
tionary state. 

We have integrated numerically Eq. l|30[l . The usual 
FFT split technique, which exactly preserves the or- 
thonormal nature of the u's and v's, is actually not satis- 
factory because it assumes that the gap function remains 
constant in time during one time step, which breaks the 
self-consistency of the equations and leads to a violation 
of the conservation of the mean number of particles. We 
therefore used an improved splitting method detailed in 
the appendix IbI 

In all the simulations that we present in this paper, 
the trap anisotropy was e = 0.1, the chemical potential 
of the initial state of the gas was fixed to /i = 8fio;; setting 
H = h 2 kp/2m, the 2D scattering length was fixed to the 
value a2D = {h/muj) 1 / 2 = aho such that fep^D = 4; the 
rotation frequency was turned on with the following law 

ft(t) = ft sin 2 (^j for < t < t (33) 

with a ramping time t = 160a; -1 much larger than the 
oscillation period of the atoms in the trap. For t > r, the 
rotation frequency remains equal to ft. The presence of 
vortices is detected by calculating the winding number of 
the phase of the gap parameter around each plaquette of 
the grid. We also calculated the total angular momentum 
of the gas. In all the simulations, we evolved the system 
for a total time of 1000 a; -1 . 

Simulations on a small 32 x 32 grid: For such a grid 
size, the calculation time remains reasonable so that we 



varied the final rotation frequency in steps of 0.1a;. For 
final rotation frequencies ft < 0.3a;, no vortices are found 
to enter the cloud and the cloud remains round. 

For ft = 0.4a;, the cloud remains round but a corru- 
gation of the surface of the cloud is observed to appear 
at time t ~ 240a; -1 ; the amplitude of the corrugation 
increases and two diametrically opposite vortices enter 
the cloud gently at t ~ 400a; -1 and, after a time inter- 
val of ~ 100a; -1 , settle in a stationary pair of vortices 
close to the trap center. At t ~ 700a; -1 , a second pair 
of vortices starts entering with the same mechanism; it 
then interacts with the first pair. The 4 vortices arrange 
in a stationary square at t ~ 850a; -1 . For ft = 0.5a;, 
the situation is similar: one vortex pair enters, then a 
second one, then a triplet of vortices starts entering at 
t ~ 490a; -1 ; eventually, at t > 610o> -1 the seven vortices 
arrange in a stationary regular pattern, consisting of an 
hexagon plus a vortex in the center. Selected images of 
the movie of the simulation for ft = 0.5a; are shown in 
figure E| For ft = 0.6a;, the scenario is slightly different. 
The corrugation of the surface is stronger, and a rect- 
angular pattern of 4 vortices enter at time t ~ 220a; -1 , 
shortly followed at t ~ 250a; -1 by a second rectangular 
pattern of 4 vortices. After an interaction period, 6 vor- 
tices align in the cloud in two parallel rows whereas two 
vortices are pushed away. Then a third rectangle of vor- 
tices enter. At later times, several extra vortices join the 
group; from t ~ 500a; -1 til the end of the simulation, 12 
vortices are present in the cloud, forming an almost sta- 
tionary and regular pattern. Clearly, in these scenarios, 
no global turbulence of the cloud is involved, since the 
first entering vortices are arranged in a preformed pat- 
tern obeying the parity symmetry of the Hamiltonian. 

For ft = 0.7a; and ft = 0.8a; the dynamics is very dif- 
ferent from the previous one. The shape of the cloud 
strongly elongates and deforms. Then strong turbulence 
sets in, at t ~ 160a; -1 for ft = 0.7a; (t ~ 135a; -1 for 
ft = 0.8a;), while the cloud anisotropy reduces, the den- 
sity profile becomes irregular, not only close the cloud 
boundary but also in the cloud center; one observes a 
quick entrance of disordered vortices in the cloud at time 
t ~ 190a; -1 for ft = 0.7a; (t ~ 150a; -1 for ft = 0.8a;); sev- 
eral anti- vortices reach the borders of the cloud for ft = 
0.7a; and even reach high density regions for ft = 0.8a;. 
After some evolution time, the density profile recovers a 
smooth and elliptic shape, the anti-vortices are expelled 
from the cloud and the vortex positions slowly relax to 
form a 17 (or 25 for ft = 0.8a;) vortex 'lattice' at times 
~ 500a; -1 , that remains essentially stationary til the end 
of the simulation. 

In conclusion, two distinct scenarios of vortex lattice 
formation are observed in the 32 x 32 simulations. For 
the lower rotation frequencies, a gentle entry of an or- 
dered pattern of vortices is observed. For the higher ro- 
tation frequencies, turbulence sets in and leads to the 
abrupt and disordered entrance of vortices and even 
anti- vortices, the regular and stationary vortex 'lattice' 
forming after some evolution time. Another difference 
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FIG. 6: For the numerical simulation of the 2D time depen- 
dent BCS equations on a 32 x 32 grid, density plots of the 
density of the trapped gas at selected times (in units of a; -1 ), 
for a final rotation frequency Q = 0.5a;. The trap anisotropy 
was e — 0.1 and the 2D scattering length a-zD = \/h/muj, 
and \x = 8hui in the initial state. The full spatial width of 
the simulation grid is shown in the figure. Crosses: positive 
charge vortices. Circles: negative charge vortices. From top 
to bottom and from left to right: t = 218oj -1 : a corruga- 
tion of the surface appears; t = 260o;~ 1 : a vortex pair enters; 
t — 342o> _1 : a second vortex pair enters; t = 400a; -1 : the 
vortices arrange on a square; t — 490oj -1 : a triplet of vortices 
starts entering; t — 610a> -1 : a stationary 7- vortex pattern. 



between the two scenarios is the temporal behavior of 
the density profile at the location of the vortex cores: 
whereas a dip in the density is visible from the start when 
a vortex enters the cloud with the gentle scenario, such 
a dip at the vortex location forms only after some relax- 
ation time in the turbulent scenario. 

The physical origin of the turbulent scenario is ex- 
pected to be the dynamic instability of the mode of de- 
gree n = 3 discussed in section [n] and the obtained 
movies qualitatively agree with that. More quantita- 
tively: for e = 0.1 a significant Lyapunov exponent is 



obtained for SI > 0.68a;; this is compatible with the 
fact that the numerical simulation observes turbulence 
for SI > 0.7a; only. 

What is the physical origin of the gentle scenario ? 
The observed corrugation at the surface suggests that it 
is driven by the instability of some surface mode local- 
ized at the surface of the cloud, which is reminiscent of 
the Landau mechanism. As a test of this idea, we have 
performed a numerical calculation of the stationary BCS 
state in a rotating frame, by the above mentioned iter- 
ative scheme: as shown in figure {7\ giving the angular 
momentum of the stationary BCS solution as a function 
of the rotation frequency, for e = 0.1, the branch with 
no vortex is followed up to SI — 0.3a;: for larger values of 
SI, the algorithm jumps to a configuration with vortices. 
This suggests that the vortex free BCS state is indeed 
not a local minimum of energy for SI > 0.3a;. What is 
then puzzling at this stage is that a harmonic stirrer can 
excite only the quadrupolar modes, whereas the negative 
energy surface mode initiating the Landau mechanism 
is expected to have a higher angular momentum |13|. A 
possible solution to this paradox was obtained by running 
the time dependent simulation with e = 0, for SI = 0.5a;: 
in this case, the harmonic trap, being isotropic, can not 
stir the gas and the stirring is due only to the fixed pe- 
riodic boundary conditions in the rotating frame; still, 
vortices were found to enter the cloud. This shows that 
the quantization box in our 32 x 32 simulations is small 
enough to activate the Landau mechanism. 

Simulations on larger grids: To get rid of the previously 
mentioned finite quantization box effects, we performed 
simulations on larger grids, 48 x 48 and 64 x 64. For the 
48 x 48 grids, we investigated the rotation frequencies 
from 0.4a; to 0.8a; in steps of O.Icj. No vortex is found 
to enter the cloud for SI < 0.5a;. For SI — 0.6a;, vortices 
enter according to the gentle scenario; the first vortices 
(in the form of a pair) enter however at a considerable 
later time, t ~ 400a; -1 , than with the 32 x 32 grid simu- 
lation. For SI > 0.7a; the vortices enter according to the 
turbulent scenario. The turbulent scenario on the 48 x 48 
grid is similar to the one observed on the 32 x 32 grid, 
except for temporal shifts: e.g. on the 48 x 48 grid, the 
turbulent period starts ~ 240a; -1 later for SI = 0.7a; and 
~ 65a; -1 later for SI = 0.8a;. 

For the 64 x 64 grids, the CPU time for a single realiza- 
tion exceeds one month on a bi-processor AMD Opteron 
workstation, so that we have considered only two val- 
ues of the rotation frequency. For SI = 0.6a;, no entry 
of vortices is observed. This confirms that the observa- 
tion of the gentle scenario, at least up to the maximal 
evolution time (1000a; -1 ) considered here, is an artifact 
of the quantization box. For SI = 0.8a;, vortices enter 
according to the turbulent scenario. The timing is now 
quantitatively the same at the 48 x 48 simulation: in both 
simulations, the vortices enter the cloud at t ~ 205a; -1 
and crystallize in a quasi stationary pattern at times 
~ 500 — 550a; -1 . Selected images of the movie of the 
64 x 64 simulation for SI = 0.8a; are shown in figure|S] 
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FIG. 7: In a steady state solution of the 2D BCS theory, 
on a 32 x 32 grid, angular momentum per particle in the 
gas, in units of h, as a function of the rotation frequency, 
for e = 0.1 and /i = 8hoj, a^o = (h/mcu) 1 ^ 2 . Black disks: 
numerical result from an iterative algorithm (no vortex at 
the left part of the jump, vortices present at the right part of 
the jump). Solid line: hydrodynamic prediction (no vortex). 
The upper graph corresponds to the same data as the lower 
graph, but for a different scale of the vertical axis: it shows 
the good agreement of the hydrodynamic prediction with the 
full numerics. In the numerics, the total number of particles 
is ~ 72 on the left part of the jump, and reaches ~ 80 on 
the most right data point; the rotation frequency is increased 
step by step from to 0.4a;, the converged state for a given Q 
being taken as an initial guess in the iterative algorithm for 
the successive value of fi. 




FIG. 8: For the numerical simulation of the 2D time depen- 
dent BCS equations on a 64 x 64 grid, density plots of the 
density of the trapped gas at selected times (in units of a; -1 ), 
for a final rotation frequency Q = 0.8a;. The trap anisotropy 
was e = 0.1 and the 2D scattering length a,2D = \J h/muj, and 
fj, = 8hu> in the initial state. The spatial width of the simula- 
tion is truncated in the figure to have approximately the same 
width as in FigHJ] Crosses: positive charge vortices. Circles: 
negative charge vortices. From top to bottom and from left 
to right: t — 196o;~ 1 : a turbulent, elongated cloud is formed; 
t = 212a; -1 : the cloud is round again, and includes a disor- 
dered mixture of vortices and anti- vortices; t = 574a; -1 : the 
vortices crystallize in a quasi-stationary pattern; t = 998a; -1 : 
slow and small shifts of some vortex positions have taken place 
with respect to the previous density plot. 



To allow for a quantitative comparison between the 
simulations for the three grid sizes, we have plotted in 
figure OH the total angular momentum of the gas as a 
function of time, for (a) fl — 0.6a; and (b) fl = 0.8a;. 
We have also given the (vortex free) hydrodynamic pre- 
diction; remarkably, this shows that the simulations give 
results close to the hydrodynamic one as long as no vor- 
tex enters the cloud, see in particular the 64 x 64 results 
for £1 = 0.6a;. To briefly address the experimental ob- 
servability of the vortex pattern, we also show in figure 
1101 a cut of the particle density (directly measurable in an 
experiment) and of the gap parameter (not directly ac- 
cessible experimentally) for the 64 x 64 simulation with 
Q = 0.8a; at a time when the vortex lattice is crystal- 
lized, this in parallel to an isocontour of the magnitude 
of the gap parameter: vortices embedded in high den- 
sity regions result in dips in the density profile, with a 
contrast on the order here of 30%. 



IV. CONCLUSION 

We have investigated a relevant problem for the present 
experiments on two-spin component interacting Fermi 
gases, the possibility to form a vortex lattice by slow 
ramping of the rotation frequency of the harmonic trap 
containing the particles. The observation of such a vor- 
tex lattice in steady state would be a very convincing 
evidence of superfluidity |3^ |. 

For a 2D model based on the BCS theory, and for the 
3D unitary quantum gas, we predict analytically, with 
the superfluid hydrodynamic equations, that the gas ex- 
periences a dynamic instability when the final rotation 
frequency is above some minimal value 0„ that we have 
calculated. This dynamic instability is very similar to the 
one discovered for a rotating Bose-Einstein condensate of 
bosonic atoms, where it was shown to lead to the vortex 
lattice formation. 

To see if this dynamic instability leads to the forma- 
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FIG. 9: Angular momentum per particle in the gas, in units 
of h, as a function of time, for a final rotation frequency (a) 
Q, — O.Guj and (b) Q — 0.8u>. Curves (from top to bottom 
in (a)): numerical simulations of the 2D time dependent BCS 
equations for the grids 32 x 32, 48 x 48 and 64 x 64, and for the 
time dependent superfluid hydrodynamic theory of section ITT1 




FIG. 10: At time t = 574w _1 of the 64 x 64 numerical simula- 
tion for f2 = 0.8ui. Left panel: isocontours of the modulus of 
the gap parameter, showing the presence of a vortex lattice; 
the x and y coordinates run from — 10ah o to +10ciho in the 
simulation but this left panel figure is truncated to a position 
interval approximately — 7aho to +7aho- Right panel: on the 
line y — — 0.627a ho , x dependence of the density p (solid line, 
in units of a^ 2 ) and of the modulus of the gap parameter 
(dashed line, in units of huj). The gap parameter was mul- 
tiplied by 2/3 for clarity. A Fourier interpolation technique 
was used in the right panel to map the 64 x 64 simulation grid 
onto a 128x128 grid. 



tion of vortices also in the case of the Fermi gases, we 
have solved numerically the full 2D time-dependent BCS 
equations, for a trap anisotropy e = 0.1. For a final 
rotation frequency f2 above the predicted Q u , we see tur- 
bulence and the subsequent fast entry of vortices. We 
conclude that the dynamic instability can indeed result 
in a vortex lattice formation. The apparent irreversibil- 
ity and energy dissipation that this seems to imply may 
be surprising at first sight, since the equations of motion 
that we integrated are purely conservative. The clue is 
probably the same as in the bosonic counterpart of these 
simulations : the spatial noise produced in the turbu- 
lent phase populates many eigenmodes (including collec- 
tive modes) of the system, and the subsequent non-linear 
evolution leads to effective thermalization of the modes. 
For fl < f2„ but for f2 larger than what we estimated 



to be the Landau rotation frequency (above which the 
vortex free superfluid is no longer a local minimum of 
energy in the rotating frame), we also see the formation 
of a vortex lattice in the simulations on the small 32 x 32 
grids, but with a gentle mechanism not involving turbu- 
lence and leading to the entrance of a pre-formed regular 
vortex pattern from the surface of the cloud. But we 
also performed simulations on larger grids: on a 64 x 64 
grid, the gentle mechanism disappears; it is therefore an 
artifact of the periodic boundary conditions that rotate 
in the lab frame and provide an artificial stirring effect. 
In a real experiment, however, we expect such a gentle 
mechanism to occur for a gas initially at finite tempera- 
ture, when the normal component of the gas is set into 
rotation by the stirrer. 
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APPENDIX A: SIMPLE DERIVATION OF THE 
HYDRODYNAMIC EQUATIONS FROM BCS 
THEORY 

We show here that the time dependent hydrodynamic 
equations Eq. 112fl and Eq. I|13l) can be formally derived for 
a vortex free gas from the time dependent BCS equations 
by using the lowest order semi-classical approximation 
and an adiabatic approximation for the resulting time 
dependent equations. As in the remaining part of the 
paper, we consider here the regime where the chemical 
potential is positive and larger than the binding energy 

The general validity condition of a semi-classical ap- 
proximation is that the coherence length of the gas 
should be much smaller than the typical length scales 
of variation of the applied potentials. Two coherence 
lengths appear for a zero temperature BCS Fermi gas: 
the inverse Fermi wave- vector, kp , associated to the 
correlation function (i/;|(r)'0|(r')), and the pair size, 
fees ~ h 2 kp/m\ A|, associated to the correlation function 
(■0l(r)?/>x(r'))- A first typical length scale of variation of 
the matrix elements in Ea. (|30|l comes from the position 
dependence of |A|: in the absence of rotation, we assume 
that this is the Thomas-Fermi radius Rtf of the gas, de- 



fined as h 2 kp/2m 



2 R 2 ?F /2. This assumes that the 



scale of variation of the modulus of the gap is the same 
as the one of the density; the adiabatic approximation to 
come will result in a |A| related to the density by Eq.Q, 
which justifies the assumption. Necessary validity condi- 
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tions of a semi-classical approximation are then: 

kp 1 ,lBCS<Rj:F. (Al) 

In the BCS regime regime, ftp 1 < Ibcs'- for an isotropic 
harmonic trap, one then finds that the condition Ea. l|Al|) 
is equivalent to 



|A| > fuv, 



(A2) 



where uj is the atomic oscillation frequency. The pres- 
ence of vortices introduces an extra length scale in the 
variation of |A|, on the order of ZbcSj which invalidates 
the semi-classical approximation. 

In the rotating case, however, this is not the whole 
story, as the phase of A can also become position de- 
pendent. As we shall see, the phase of A in this pa- 
per may vary as mtuxy/h: when this quantity varies by 
~ 2ir, A changes completely; this introduces a length 
scale ~ 2nh/(mujR^F) ~ 1/kp, making a semi-classical 
approximation hopeless. We eliminate this problem by 
performing a gauge transform of the it's and u's: 

fi.(r,f) = u s (r,t)e- lS ^/ h (A3) 
v s (r,t) ee v s (r,t)e +tS ^/ h (A4) 

where the phase is defined in Eci. (|10fl . The time depen- 
dent BCS equations are modified as follows: 



ihd t 



h |A| 
IAI -U 



EE L 



(A5) 



where the gauge transformed Hamiltonian is 

ho = e~ lS l h hoe +lS/h + d t S. (A6) 

Let us review relevant observables in the gauge trans- 
formed representation. First the gap equation is modified 

as 



|A| = -go^2u 8 v* a 



Then the mean total density reads 
P = 2^u s u*. 



(A7) 



(A8) 



Last, we introduce the total matter current j(r,i), that 
obeys by definition 



d t p + div j = 0. 



(A9) 



In the rotating frame, in a many-body state invariant by 
exchange of the spin states f and J,, it is very generally 
given by 

j = — (V^lgradiM — c.c.) — pSi x r. (A10) 
im V 1 / 

Within BCS theory, this gives 
j = /)(v-flxr)| - [v* grad v s — v s grad v*] , 

TO Z — ' 

(AH) 



where the velocity field v is defined as grad S/m. Note 
that the continuity equation Eq. (|A9() remains true for the 
BCS theory j^. 

To calculate the two key quantities Ea. (|A8|l and 
Ea. (|All|l . it is sufficient to know the following one-body 
density operator for a fictitious particle of spin 1/2, 



G= ( ^TT <rn \ = ( \u*s)(u s \ \u s )(v s \ 

°"IT a u J ^ V I s *) l^)(^«l 



(A12) 



To prepare for the semi-classical approximation we in- 
troduce the Wigner representation of a [38|: 

W(r,p,*)=Wigner{a}EE J A^ 3 ( r -x/2\a\v+x/2)e i ^/ H 

(A13) 

where d is the dimension of space. The key observables 
have then the exact expressions: 

p(r,t) = 2 J d d pW Ll (r,p,t) (A14) 

|A|(r,t) = -go J d d pW n (r,p,t) (A15) 
j(r,t) = p(v - n x r) 

~ f d d ppW u (T,p,t). (A16) 

The semi-classical expansion then consists e.g. in 

Wigner{l/(f)cj} = [V(r) + ^d T V ■ d p + . . ]W(r, p, t). 

(A17) 

The successive terms we called zeroth order, first order, 
etc, in the semi-classical approximation. 

We write the equations of motion Eq. I|A5|I up to zeroth 
order in the semi-classical approximation: 

ihd t W(r,p,t)\W = [£ (r,p,t),W(r,p,f)] (A18) 
where the matrix Lo is equal to 

_ d-^n{v,t) |A|(r,t) \ 



Lo(r,p,t)= 2 



|A|(r,<) -i^ + Meff (r,t) 



(A19) 

We have introduced the position and time dependent 
function, 

/i c ff(r, t) = /x- U(r,t) - ^mv 2 +mv (SI x r) - d t S(r, t), 

(A20) 

that may be called effective chemical potential for reasons 
that will become clear later. 

At time t = 0, the gas is at zero temperature. By 
introducing the spectral decomposition of Lit = 0) one 
can then check that 



a(t = 0) = 6[L(t = 0)] 



(A21) 



where 9(x) is the Heaviside function. Since Lo(t — 0) is 
the classical limit of the operator Lit = 0), the leading 
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order semi-classical approximation for the corresponding 
Wigner function is, in a standard way, given by 



W(r,p,f = 0) 



1 



{2irh) 



3 0[£o(r,p,t = O)] (A22) 



that is each two by two matrix W is proportional to a 
pure state IV") (^| with 



|V(r,p,i = 0)) 



_ fU (v,p) 
V (r,p) 



(A23) 



where (Uo, Vq) is the eigenvector of £ (r, p, t = 0) of posi- 
tive energy and normalized to unity. At time i, according 
to the zeroth order evolution Ea. (|A18|l . each two by two 
matrix W remains a pure state, with components U and 
V solving 



ihdt 



t/(r,p,t) 
V{r,p,t) 



= L Q (r,p,t) 



(U{v,p,t) 
\V(r,p,t) 



(A24) 



We then introduce the so-called adiabatic approxima- 
tion: the vector (U, V) , being initially an eigenstate of 
Lo(r,p,t = 0), will be an instantaneous eigenvector of 
Lq{t,p,£) at all later times t. This approximation holds 
under the adiabaticity condition 39], detailed below, re- 
quiring that the energy difference between the two eigen- 
values of L (r,p,t) (divided by H) be large enough. As 
this energy difference can be as small as the gap param- 
eter, this will impose a minimal value to the gap, as we 
shall discuss later. In this adiabatic approximation, one 
can take 



W(r,p,t) 



1 



6[L (r,p,t)] 



1 



(2irh) d L uv n (2nK) d 

(A25) 

where | + (r,p,i)), of real components (J7i„ s t, Mnst), is 
the instantaneous eigenvector with positive eigenvalue of 
the matrix Lo defined in Ea. (|A19f) . Its components are 
simply the amplitudes on the plane wave exp(ip • r/h) 
of the BCS mode functions of a spatially uniform BCS 
gas of chemical potential /i e g and of gap parameter 
|A(r, t)|. Using Ea. HA14|) and Ea. (|A15fl with the approx- 
imate Wigner distribution Ea. (|A25(l . one further finds 
that this fictitious spatially uniform BCS gas is at equilib- 
rium at zero temperature so that the expressions Eq.J2J 
and Eq.Q) may be used. In particular, Eq.JSJ) gives 



ju e ff(r,i) = no[p(r,i)] 



(A26) 



which leads, together with Ea. (|A20|l . to one of the time 
dependent hydrodynamic equations, the Euler-type one 
Ea. (jT3")l . Also, ?7inst and Vinst are even functions of p, 
so that the integral in the right hand side of Eq. l|A16|) 
vanishes and Ea. i|A9() reduces to the hydrodynamic con- 
tinuity equation Eq. l|12|) . Under the adiabatic approxi- 
mation, the superfluid hydrodynamic equations are thus 
derived. 

We now discuss the validity of the adiabatic approx- 
imation. Without this approximation, the two by two 



matrix W has non-zero off-diagonal matrix elements 
(+| W| — ) where |— ) is the instantaneous eigenvector of 
Ea. (|A19|l with a negative eigenvalue, that can be written 
(Mnst, — CAnst)- Writing from Eo. (|A18() the equation of 
motion for (+|W|— ), one indeed finds a coupling to the 
diagonal element (+|W|+) due to the non infinite ramp- 
ing time of the rotation. This coupling can be calculated 
using the off-diagonal Hellman-Feynman theorem for real 
eigenvectors, and corresponds to a Rabi frequency 



1 



r^tir 



-\dt 



1 



(d t L ) 



(A27) 



where e± is the eigenenergy of |±) for the matrix Lq 



e± =± (p 2 /(2m) -/icff) 2 + |A| 5 



1/2 



(A28) 



But this is not the whole story, as we have neglected 
the so-called motional couplings, that can also destroy 
adiabaticity. These motional couplings are due to the 
fact that |+) and |— ) depends on r, p and that terms 
involving d p W and d r W will appear in the equation for 
W beyond the zeroth-order semi-classical approximation. 
Such non-adiabatic effects are well known for the motion 
of a spin 1/2 particle in a static but spatially inhomo- 
geneous magnetic field. In our problem, the first order 
term of the semi-classical expansion is actually simple to 
write: 

d t W\ {1) = - [d r L ■ d p W - d p L ■ d r W + h.c] . (A29) 

The matrix L corresponds to the classical limit of L(t): 

L(r, p, t) = L (r, P, *) + P • (v - Vt x r) I, (A30) 

where I is the 2x2 identity matrix. In the resulting 
equation of evolution of (+|W| — ), taking (+|W|+) = 
l/(2irh) d and (— \W\— ) = 0, a motional Rabi coupling to 
(+| W]+) now appears: 

^motion = -~d p [p ■ (v - n x r)] ■ (-\d T \+) 

+a r [ P -(v-nxr)]-<-|a p |+)(A3i) 

Expressions similar to the one for (— |Qt|+) can be derived 
with the off-diagonal Hellman-Feynman theorem. 

We now calculate the total Rabi frequency ftot = 
I'time + ^motion at the local Fermi surface, that is for a 
value of the momentum such that p 2 /2m — (x e fi(r,t). 
This is indeed at the Fermi surface that we expect the 
adiabaticity condition to be most stringent, as the energy 
difference e+ — e_ takes there its minimal value, equal to 
twice the gap |A(r,i)|. Then Ui nst = V5 nst = 1/V2 and 
the expressions resulting from the Hellman-Feynman the- 
orem are very simple: 



■\d x 



9\(jieS -p 2 /2m) 
2IAI 



(A32) 
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where A stands for t or for an arbitrary component of the 
vectors r or p. We then get the condition for adiabaticity: 



\Vtot\ 



2IAI 



Dt 



S 



<2|A|/fi, (A33) 



where D/Dt = d t + (v — tt x r) ■ d r . 

A fully explicit expression for i/ to t can be obtained us- 
ing the hydrodynamic equations and taking the limit of 
a very long ramping time of the rotation, as is the case 
in our simulations, so that the hydrodynamic variables 
are close to a steady state and S ~ mui/3(t)xy. Using 
Ea. (|A26|l and the continuity equation Ea. (|12fl . one gets 
Dfics/ 'Dt = — p[i' Q [p]divv ~ so that one is left with 



1 



-Vtot 



/3{t)u>P x Py 



(A34) 



The constraint 
dition in 2D: 



"tot /2 1 -C 2\A\/h then results in the con- 



hu < 4£V|/3(t)|, 



(A35) 



evolving the (u' s ,v' s ) with the r-dependent part of two 
by two matrix of Ea. l3U|) during dt, for a fixed value of 
A(r,i) = -jo^ s !is(r)»5*(r). This exactly preserves the 
unitary of the full evolution, but the fact that a fixed 
value of A is taken during the second step of the evolu- 
tion breaks the self-consistency between A and u s , v s so 
that the total number of particles, N — 2 ^2 s (v s \v s ), is 
conserved to first order in dt but not to all orders in dt. 
Numerically, for the time steps dt leading to a reason- 
able CPU time, one then observes strong deviations of 
this total number from its initial value. Note that such 
a problem does not arise for the time dependent Gross- 
Pitaevskii equation for bosons, for which conservation of 
unitary and number of particles is one and a same thing. 

This problem for the BCS equations can be fixed by 
restoring the self-consistency for the evolution during dt 
associated to the r-dependent part of the equation of 
evolution. That is one solves during dt: 



ihdt 



U(r) 
A*(r 



-M A(r,f) 
t) h-U{t) 



(Bl) 



where Eq is the dimer binding energy. This condition 
is satisfied in our simulations as f3 is at most ~ 0.64 
(for = 0.8a;) and we took a 2 D = (fr/muj) 1 / 2 , \i = 8hio 
resulting in Eq ~ l.ihus and A ~ i.7fko. Note that it is in 
general more stringent than the usual condition Eq. i|A2|) 
but for the particular parameters of our simulations, it 
turns out to be roughly equivalent. 



not for a fixed A but with the time dependent A given by 
the self-consistency condition Eq.J2SJ. As a consequence, 
Ea. (IBl|) written for all modes s is a set of non- linearly 
coupled time dependent equations. Fortunately, they are 
purely local in r, so that they can be solved analytically. 
One finds that A(r,t) varies as e ~ lX ( r ) t j where 



APPENDIX B: A SPLITTING TECHNIQUE 
CONSERVING THE MEAN NUMBER OF 
PARTICLES 

The standard splitting technique approximates the 
evolution due to Ea. (|30|l during a small time step dt 
by first evolving the (u s ,v s ) into (u' s ,v' s ) with the ki- 
netic energy and rotational energy during dt, and then 



H\(r) = 2[U(r) - fi] 



g J2[\<r,t)\ 2 -\u s (r,t)f 



(B2) 

can be checked to be time independent for the local 
evolution Ea. ljBlfl . Then the system Ea. i|Bl|) is trans- 
formed into one with time independent coefficients (so 
readily integrable) by performing a time dependent gauge 
transform, u s (r,t) = U s {T,t)e^ lX ^ t ^ 2 and v s (r,t) = 
V s {v,t)e +lX ^ t / 2 . 
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